# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
# install.packages(c("tidyverse", "estimatr", "stargazer"))

library(tidyverse)
library(estimatr)
library(stargazer)

# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")
source("opeds_source.R")

dv_df <- dv_df %>% filter(issue != "Climate")
mturk_opeds <- filter(mturk_opeds, Z != "climate")

mturk_opeds <-
  mturk_opeds %>%
  mutate(Z_match = factor(
    Z,
    levels = c("control", "amtrak", "paul", "veterans", "wallstreet"),
    labels = c("control", "amtrak", "flat", "vets", "wall")
  ))

elite_opeds <-
  elite_opeds %>%
  mutate(Z_match = factor(
    Z,
    levels = c("control", "amtrak", "paul", "veterans", "wallstreet"),
    labels = c("control", "amtrak", "flat", "vets", "wall")
  ))

table(mturk_opeds$Z)
table(elite_opeds$Z)
table(mturk_opeds$Z_match)
table(elite_opeds$Z_match)

stacked_opeds <-
  bind_rows(mturk = mturk_opeds, elite = elite_opeds, .id = "experimental_sample") %>%
  mutate(experimental_sample = factor(experimental_sample, levels = c("mturk", "elite")))

# Interaction Term --------------------------------------------------------

fit_1 <- lm(dv_amtrak_scale_w1 ~ Z_match*experimental_sample, data = filter(stacked_opeds, responded_w1))
fit_2 <- lm(dv_flat_scale_w1 ~ Z_match*experimental_sample, data = filter(stacked_opeds, responded_w1))
fit_3 <- lm(dv_vets_scale_w1 ~ Z_match*experimental_sample, data = filter(stacked_opeds, responded_w1))
fit_4 <- lm(dv_wall_scale_w1 ~ Z_match*experimental_sample, data = filter(stacked_opeds, responded_w1))


stargazer(
  fit_1,
  fit_2,
  fit_3,
  fit_4,
  se = starprep(fit_1, fit_2, fit_3, fit_4),
  p = starprep(fit_1, fit_2, fit_3, fit_4, stat = "p.value"),
  style = "apsr",
  column.sep.width = "0pt",
  digits = 3,
  omit.stat = c("adj.rsq", "f", "ser"),
  dep.var.labels = c("Amtrak", "Flat Tax", "Veterans", "Wall Street"),
  covariate.labels = c(
    "Op-ed: Amtrak",
    "Op-ed: Flat Tax",
    "Op-ed: Veterans",
    "Op-ed: Wall Street",
    "Elite Experiment",
    "Elite X Amtrak",
    "Elite X Flat Tax",
    "Elite X Veterans",
    "Elite X Wall Street",
    "Constant (Constant)"
  ),
  model.numbers = FALSE,
  notes = c(
    "Models estimated via OLS. Robust standard errors in parentheses.",
    "Dependent variables are constructed via factor analysis and have unit variance."
  ),
  notes.append = TRUE,
  font.size = "small",
  float = TRUE,
  align = TRUE,
  label = "tab:interaction_scale",
  title = "Comparison of Treatment Effects on Composite Scale Dependent Variables"
)


# Make Scatterplot --------------------------------------------------------

df <- NULL

for(i in 1:nrow(dv_df)){
  formula_w1 <- formula(paste0(dv_df$dvs[i], "_w1 ~ Z_match"))
  formula_w2 <- formula(paste0(dv_df$dvs[i], "_w2 ~ Z_match"))
  fit_1 <- tidy(lm_robust(formula_w1, data = filter(mturk_opeds, responded_w1 == TRUE)))
  fit_2 <- tidy(lm_robust(formula_w1, data = filter(elite_opeds, responded_w1 == TRUE)))
  fit_3 <- tidy(lm_robust(formula_w2, data = filter(mturk_opeds, responded_w2 == TRUE)))
  fit_4 <- tidy(lm_robust(formula_w2, data = filter(elite_opeds, responded_w2 == TRUE)))
  
  df_w1 <- 
    data.frame(
      mturk_est = fit_1[,"estimate"],
      mturk_se = fit_1[,"std.error"],
      mturk_conf.low = fit_1[,"conf.low"],
      mturk_conf.high = fit_1[,"conf.high"],
      elite_est = fit_2[,"estimate"],
      elite_se = fit_2[,"std.error"],
      elite_conf.low = fit_2[,"conf.low"],
      elite_conf.high = fit_2[,"conf.high"],
      dv = dv_df$dvs[i],
      dv_label = dv_df$dv_labels[i],
      dv_cat = dv_df$dv_cat[i],
      wave = "Wave 1",
      treatment = fit_1[,"term"]
    )
  
  df_w2 <- 
    data.frame(
      mturk_est = fit_3[,"estimate"],
      mturk_se = fit_3[,"std.error"],
      mturk_conf.low = fit_3[,"conf.low"],
      mturk_conf.high = fit_3[,"conf.high"],
      elite_est = fit_4[,"estimate"],
      elite_se = fit_4[,"std.error"],
      elite_conf.low = fit_4[,"conf.low"],
      elite_conf.high = fit_4[,"conf.high"],
      dv = dv_df$dvs[i],
      dv_label = dv_df$dv_labels[i],
      dv_cat = dv_df$dv_cat[i],
      wave = "Wave 2",
      treatment = fit_3[,"term"]
    )
  
  df <- bind_rows(df, df_w1, df_w2)
}

df <- df %>%
  mutate(oped_dv_match = ifelse(gsub("Z_match", "", treatment) == as.character(dv_cat), "Target Issue", "Non-Target Issue")) %>%
  filter(treatment != "(Intercept)")


# These figures not in paper or appendix but were in presentations --------

g <-
  df %>% filter(wave == "Wave 1") %>%
  ggplot(aes(
    x = mturk_est,
    y = elite_est,
    color = oped_dv_match,
    shape = oped_dv_match
  )) +
  geom_point(aes(alpha = oped_dv_match)) +
  geom_vline(xintercept = 0,
             linetype = "dashed",
             alpha = .5) +
  geom_hline(yintercept = 0,
             linetype = "dashed",
             alpha = .5) +
  geom_segment(aes(
    x = mturk_est,
    xend = mturk_est,
    y = elite_conf.high,
    yend = elite_conf.low
  ),
  alpha = .2) +
  geom_segment(aes(
    x = mturk_conf.high,
    xend = mturk_conf.low,
    y = elite_est,
    yend = elite_est
  ),
  alpha = .2) +
  theme_bw() +
  coord_cartesian(xlim = c(-.3, .8), ylim = c(-.3, .8)) +
  scale_alpha_manual(values = c(.3, 1)) +
  scale_color_manual(values = c("blue", "red")) +
  xlab("Mechanical Turk Sample Standardized Treatment Effect Estimates") +
  ylab("Elite Sample Standardized Treatment Effect Estimates") +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    strip.background = element_blank()
  )
g



g <-
  ggplot(df,
         aes(
           x = mturk_est,
           y = elite_est,
           color = oped_dv_match,
           shape = oped_dv_match
         )) +
  geom_point(aes(alpha = oped_dv_match)) +
  geom_vline(xintercept = 0,
             linetype = "dashed",
             alpha = .5) +
  geom_hline(yintercept = 0,
             linetype = "dashed",
             alpha = .5) +
  geom_segment(aes(
    x = mturk_est,
    xend = mturk_est,
    y = elite_conf.high,
    yend = elite_conf.low
  ),
  alpha = .2) +
  geom_segment(aes(
    x = mturk_conf.high,
    xend = mturk_conf.low,
    y = elite_est,
    yend = elite_est
  ),
  alpha = .2) +
  theme_bw() +
  coord_cartesian(xlim = c(-.3, .8), ylim = c(-.3, .8)) +
  scale_alpha_manual(values = c(.3, 1)) +
  scale_color_manual(values = c("blue", "red")) +
  facet_wrap( ~ wave) +
  xlab("Mechanical Turk Sample Standardized Treatment Effect Estimates") +
  ylab("Elite Sample Standardized Treatment Effect Estimates") +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    strip.background = element_blank()
  )
g

# Correlations
df %>%
  group_by(oped_dv_match, wave) %>%
  summarize(cor(mturk_est, elite_est),
            p = cor.test(mturk_est, elite_est)$p.value)


test <- with(df, cor.test(mturk_est, elite_est))
test$p.value

